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Abstract 

A new model to describe fractal growth is discussed which includes effects due to 

long-range coupling between displacements u. The model is based on the biharmonic 

O 

Q \ equation V^m = in two-dimensional isotropic defect-free media as follows from the 

O I Kuramoto-Sivashinsky equation for pattern formation -or, alternatively, from the theory 

of elasticity. As a difference with Laplacian and Poisson growth models, in the new 
t~-~. ■ model the Laplacian of u is neither zero or proportional to u. Its discretization allows 

O ■ to reproduce a transition from dense to multibranched growth at a point in which the 
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growth velocity exhibits a minimum similarly to what occurs within Poisson growth in 
planar geometry. Furthermore, in circular geometry the transition point is estimated for 
the simplest case from the relation r^ ^ L/e^'"^ such that the trajectories become stable 



'^ ' at the growing surfaces in a continuous limit. Hence, within the biharmonic growth 



model, this transition depends only on the system size L and occurs approximately at a 
distance 60% far from a central seed particle. The influence of biharmonic patterns on 
the growth probability for each lattice site is also analysed. 
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1 Introduction 

The study of pattern formation in physically important fields by using simple lattice 
models under certain boundary conditions is of great interest because they describe many 
phenomena as occur in nature. At present, numerical simulation studies include, e.g., the 
dielectric breakdown model (DBM), i.e., solving the Laplacian equation in the medium 
surrounding the growing aggregates [|l], @, |^ , and the Poisson growth model to the pattern 
formation in screened electrostatic fields 0]. Besides, the diffusion-limited aggregation 
(DLA) model also plays a crucial role in illustrating fractal growth [0, ^, |^. So far, many 
efforts have been invested in such stochastic models that share the common feature of 
being essentially second-order differential equations. In these models iterative procedures 
are carried out around four mesh points of a lattice to then phenomenologically relate 
the global influence of a growing pattern to the growth probability for each lattice site 
under different power-law forms (see, e.g., [|[). 

Nevertheless, there is also another important class of physical problems leading to 
partial differential equations that can also be linear in the order parameter but, as a 
peculiarity, they are of higher order and directly concern the issue of pattern formation. 
An example of this is the (time-averaged) Kuramoto-Sivashinsky equation |^, |TD[: 

^ = ^,^^u + \W^W^u + ^^{Wuf , (1) 

which models pattern formation in different physical contexts, such as chemical reaction- 
diffusion systems [|l5, [l^ and a cellular gas flame in the presence of external stabilizing 



factors [13|. In the above z/, A and /i are normalized coefficients and u a displacement 
at time t. 

If we consider the simplest case, i.e. assume static solutions and keep linear terms 
only {i.e., fi = 0), then Eq.(|I|) reduces to the sum of two important terms: the Laplacian 



equation plus the Biharmonic equation 

VS = V\V^u) = . (2) 

Since the former term [i.e. A = 0), when discretized, is known to display fractal struc- 
tures on a squared lattice [jl], 0, then the latter term {i.e., Eq.(^)) might also in fact 
select a new class of fractal patterns. 

In addition to the above physical relevance, the biharmonic part of Eq.(0) can de- 
scribe the deflection of a thin plate subjected to uniform loading over its surface with 
fixed edges §, Q, the steady slow two-dimensional (2D) motion of a viscous fiuid [l^ 



or the vibration modes in the acoustic of drums [0. Besides this a higher-order dif- 
ferential equation containing the biharmonic term also appears in the study of kinetic 
growth with surface relaxation (see, e.g., fl^, [18|, |T9|, pO|). 



As an important difference with respect to second-order, to solve numerically higher- 
order differential equations, such as the one at hand, requires values for the order pa- 
rameters at either their first or second normal derivatives at each boundary point and 
beyond |2T|, ^. Hence, an analysis of effects due to long-range (many-body) coupling of 



lattice sites, including those mesh points at the lattice boundaries, on the formation of 
connected patterns is by no means trivial and as such it is a completely open problem. 
This is the subject of this paper. 

In the following we shall discuss what we believe to be the first attempt to in- 
clude higher nearest neighbour (displacements) shells in numerical simulations of fractal 
growth in isotropic defect-free media of arbitrary elastic constants. We achieve this by 
(analytically and numerically) studying Eq.(0) under the condition that the Laplacian 
of M is not constant {i.e., Laplace's case) or proportional to u {i.e., Poisson's case) along 
a growing surface. Since this biharmonic operator can also be generated in the theory 
of elasticity ||14|, ^, ^ , then u = r — r' may become the displacement vector and the 



deformation of a body may be caused by an applied force which appears in the solution 
of Eq.(0) through some boundary conditions (Eq.(^) below). 

2 Numerical simulation 

In this work, the numerical simulation of Eq.(^ is carried out on a lattice with either 
planar {i.e., growth in a channel) or circular {i.e., radial growth) geometries. Within 
the planar geometry we use a 100 x 200 lattice, having periodic boundary conditions 
along the x-direction, and set the values u° = 1 and m* = for the upper ?/-boundary 
and the surface aggregate {y = £), respectively. Initially, seed particles are placed on a 
line at displacement u{i, 1). For circular geometry, we use lattice sites enclosed within a 
circunference of radius r = 100 and locate only one seed at the center such that u° and 
■u* are also unity and zero at the outer (circular) boundary and at the inner (growing) 
aggregate, respectively. 

The inclusion of (second and higher-) nearest neighbour bond shells in our numerical 
simulations follows the discrete form of Eq.(^ on the {i,j) lattice site, which yields an 
expression involving values of u at thirteen mesh points written as pi 



Ui-2,j + 2ui-ij^i - Sui^ij + 2ui^ij+i + Uij-2 - Suij^i + 20uij (3) 

It is important to mention that numerical simulations using this equation are somehow 
more involved than for Laplace [0] or Poisson Q] growth because of long-range coupling. 
However, the accuracy of the solution for Eq.(^ can similarly be improved by looking 
at the convergence of the iterative solution using the Gauss-Seidel method |^ . 



Equation (|^) requires modificiation when applied at mesh points that are adjacent 
to a boundary, since one (at the edge) or two (near the corners) of the values needed are 
at sites outside the lattice. This modification is made by introducing a fictitious mesh 



point at {i, L) outside the planar lattice in the y-direction, where the value of u is there 
given by the derivative boundary condition along one edge boundary: 

, du . , , 

Ui^L+l = Ui^L-l + 2h— \ix ■ (4) 

Here h {= h^ = hy) is the mesh size which we set equal to unity for simplicity. For 

planar geometry, we evaluate Eq.(^ approximating 

du , 3(m° — u^) ,^. 

dy L 

for all 2-columns. This expression is obtained after taking the limit i/L << 1 in the 

solution of a one-dimensional biharmonic equation. Also for simplicity, the values of 

u{i, 0) for the points below the seed will be set to zero. For circular geometry we treat 

the boundary conditions similarly to Laplacian growth |Q, i.e., all the lattice sites outside 

a circle of radius r are set equal to u°, hence, we approximate for the sake of simplicity 

du/dy |y=L~ in the absence of applied field. 

The procedure for growing our biharmonic aggregates follows standard techniques 
^. First, Eq.(ll) is solved iteratively till the solutions converge to a given accuracity 
(of the order 3 x 10^^ or smaller). Second, after adopting a growth probability law, 
the aggregates stochastically grow at one (or more) perimeter(s) under a given relation 
between u and growth probability P. 

In order to have biharmonic pattern selection we shall assume that P at the grid site 
{i,j) depends on two different phenomenological (normalized) forms, namely 

where the sum runs over all of the nearest neighbour sites to an aggregagate. Model 
I implies that Pij is proportional to the local potential, similarly to DBM |]T], ^, ^. 
Whereas in model //, Pij is proportional to displacements around the {i,j) site. 
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3 Results and discussion 

As already mentioned we attempt to study effects of long-range coupling of displacements 
in simulations of fractal growth. To achieve this, we start analysing the numerical solu- 
tions of Eq.(0) in planar geometry. In Figs. 1(a), (b) and (c) we show three characteristic 
stages of growing for a biharmonic pattern by attaching one particle at each step and 
assuming model I for the probability P^. Figure 1(a) includes 2100 particles, (b) 4164 
particles and (c) 4355 particles added to the biharmonic cluster. In these numerical sim- 
ulations we estimate the derivative boundary condition (c./., Eq.(D) from the analytical 
solution given by Eq.(|) and fix m* = and (rescaled) 0° = V^m \y=L^ Qu° / L"^ such that 
u° = 1. Clearly, at a distance separation of about C. ~ 160, Fig. 1(c) displays features of a 
transition from dense to multibranched growth similarly to what occurs within Poisson 
growth in planar geometry [Q. This structure demostrates that long-range coupling of 
displacements drives such phenomena. 

Surprisely, such a transition also appears when we consider model II for Pij as can 
be seen in Fig. 2 for a biharmonic cluster containing 3428 particles. This means that 
the influence of ramified biharmonic patterns on the growth probability for each lattice 
site is on the type of the (far from equilibrium) pattern obtained and not on locating 
the transition. We add that several branches may also develop within our biharmonic 
growth model by attaching -simultaneously, and stochastically- more than one particle 
at each time step similarly to Louis et al.'s work 0]. We believe that for Poisson growth 
this transition is the result of including many-body contributions via screening in a sort 
of mean field approach. In our biharmonic model, long range coupling appears naturally 
as a consequence of dealing with Eq.(|^). As we shall show next for the case of planar 
growth, the transition point within both Poisson and biharmonic models appears when 
the growth velocities exhibit a minimum. 



In Fig. 3 we plot results for the grow velocity v along the y-direction for biharmonic 
growth under model I (curve A), i.e., proportional to the averaged value of P equal 
to V^M to consider explicitly nearest neighbour sites, and under model II (curve B), 
where it has been set proportional to the local displacement u. For comparison, we 
also include in this figure results for Poisson growth: i.e.., V^0 = /i^ [^, such that v is 
set proportional to the field | 0jj- — 0* | following Ref.Q]. The Poisson growth model 
incorporates screening (curve C): i.e., /i = A^ and antiscreening (curve D): i.e., fi = — A^. 

In this plot it can be seen that the transition from dense to multibranched growth 
coincides with the fact that v on the growing surface presents a minimum. This was first 
pointed out by Louis et al. f^ for the case of screening (curve C) (which also applies 
for antiscreening as seen in curve D). However, from the present findings, we can argue 
that this is also true for biharmonic growth (curves A and B) independently of how 
we have related the probability Pij to u{i,j). For Laplacian growth this phenomenon 
does not apper since the trend is to generate a single tip at faster velocity than in the 
cases of Poisson or biharmonic growth. It is worthwhile to point out that biharmonic 
patterns below the transition point, i.e., within the dense region (see Figs. 1 and 2), 
become not that denser as for screening [Q but are rather denser than for antiscreen- 
ing (not shown). Indeed, this can easily be understood from Fig. 3 since for Poisson 
growth (screening), such that fi > 0, the transition occurs at smaller velocities than for 
biharmonic growth. Henceforth, a Eden-like pattern can be generated due to screening 
1^. Above the transition multibranched fractals appear in all models considered, but 
within our biharmonic model the transition point depends on the system size only (as it 
will become clear below). We remark that both Poisson and biharmonic models lead to 
obtain growth velocities along the y-direction with parallel slopes above and below their 
respective transition points. We feel this to be a particular feature of fractal growth due 
to long-range couphng of displacements. 
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Let us focus next on effects due to long-range coupling on the fractal growth in 
circular geometry. This will allows us to make calculations of the fractal dimension 
for the biharmonic patterns by the standard box counting method |^. To generate 
and select biharmonic fractal structures we set, for simplicity, the derivative boundary 
condition required along the r-direction equal to zero. Figures 4(a) and (b) show the 
results for circular biharmonic growth using again models / (where P is set proportional 
to V^M, which corresponds to the potential in |^) and //, respectively. Figure 4(a) 
includes 2326 particles whereas Fig. 4(b) 2172 particles. It can immediately be seen that 
above a certain transition point, say re, these new plots give another demostration that 
long-range coupling is the most relevant aspect for the transition indepedently of the 
geometry adopted or of the given relation between Pij and u{i,j). 

The transition point ri can, to a good approximation, be estimated from the con- 
tinuous limit of Eq.(0) in cylindrical coordinates, such that z = for all polar angles 6, 
namely 

-2^-0^) - -(^) + 2(^) + r{^) = , (6) 

whose solution (far from the origin) is 

u = A + Br"^ + C \ogr + Dr'^ \ogr . (7) 

Clearly, two of the four coefficients A,B,C and D may be determined by the boundary 
conditions of the problem, i.e., by assuming that at r = rf. u{r() = m* = 0, whereas at 
r = L, u{L) = u° = 1. However, as discussed above, in addition to this we also have 
the condition required by simulations (c./., Eq.(^)): du/dr \r=L= 0, which leads to the 
following relation between B, C and D: 
B C/D ,, ^ 1 



This in conjuction with one possible condition of stabihty for the trajectories at the 
growing surfaces, that is obtained from du/dr |r=r^~ 0: i.e., B/D = — (logr^ + 1) and 
C = 0, gives 

L 
r^^ —- = 0.607 X L . 9 

Hence, we have that the transition point only depends on the system size L and occurs 
approximately at a distance (in lattice units) about 60% far from the seed particle as 
indicated by arrows in Figs. 4(a) and (b). Nicely, our simplest prediction is in accord 
with the numerical simulations. But if the derivative boundary condition du/dr \r=L is 
a constant different from zero, then Eq.(§) changes by a factor which, together with the 
above stability conditions, implies a shift of the transition length and a more complicated 
relation between ri and L. 

To end this section, we examine the fractal properties of the simulated biharmonic 
patterns displayed in Figs. 4(a) and (b), by simply counting the number of particles, N{r), 
inside a circle of increasing radius r (in lattice units) around a seed particle at the origin 
until we reach a distance rg at which the transition appears (c./., Eq.(|^)). We then plot 
it as a function of r in a log-log plot, as despicted in Fig.5, for models / (triangles) and 
// (squares) which we compare to Laplacian growth (circles). Over a decade, we obtain 
lines with slopes larger than unity and smaller than the space dimension. Thus, to a first 
approximation, the fractal dimension of our biharmonic patterns in the denser region 
approaches the value due to DBM (and DLA) within error bars ^. This illustrates the 
fractal nature of the biharmonic patters. Above r^, there is a change of slope due to 
the transition from dense to multibranched growth on contrary to the curve from the 
Laplacian model (circles) which continues to be linear indicating thus that the cluster 
grows dense. 



4 Concluding remarks 

In this work, we have shown that the discretization of the biharmonic Eq.(^ allows 
to reproduce a transition from dense to multibranched growth at a point in which the 
growth velocity v along the y-direction exhibits a minimum similarly to what occurs 
within Poisson growth in planar geometry. The physical basis for w, as plotted in Fig. 3, 
follows similarly to the dielectric breakdown model (model /: curve A), namely assuming 
V to be proportional to the averaged value of P equal to the potential V^m -to thus 
consider explicitly nearest neighbour sites. Whereas, model H: curve B, where v has 
been set proportional to the local displacement m, is simply a mathematical model. On 
the other hand, we also discussed that in circular geometry the transition point can 
be estimated from the relation r^ «i Lje^l'^ such that the trajectories become stable 
at the growing surfaces in a continuous limit. This is in reasonable agreement with 
present numerical simulations. Hence, we conclude that the transition from dense to 
multibranched growth within a biharmonic approach depends on the system size L only. 
The transition obtained from numerically solving the biharmonic equation might not 
be necessarily similar to the novel and intriguing Hecker transition (see e.g. PB| , Q). 



But, to this end, we notice that when assuming in the Kuramoto-Sivashinsky relation 
(HI), i^ < and A > 0, then we obtain (see also P]) 

^ = -V2v + VVv + 2vVv , (10) 

at 

in which v = Vu and Curl v = 0. This looks somewhat like the Navier-Stokes equation 

(for a potential flow with negative viscosity) which may be somehow related to the recent 

analysis in Ref . [^ of electrochemical deposition. 

Finally, it is important to emphasize again that -as a crucial difference to those 

models satisfying second order differential equations- to select fractal structures via the 

biharmonic part of Eq. (|l]) requires values for the order parameters (at either their first or 
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second normal derivatives at each boundary point) and the (second and higher-) nearest 
neighbour bond shells on the {i,j) lattice site which yields an expression involving values 
of u at thirteen mesh points. Thus, the formation of connected patterns within the 
biharmonic equation is not trivial at all. Because of this, our numerical simulations 
are more involved than for Laplace (or Poisson) fractal growth due to the long-range 
coupling between diplacements. 
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Figure captions 

• Figure 1: Three characteristic stages of growing for a biharmonic pattern by 
attaching one particle at each step and assuming model / for the probability Pj^: 
(a) 2100 (b) 4164 and (c) 4355 added particles to diplay the transition from dense 
to multibranched growth. 

• Figure 2: Final biharmonic pattern above the transition when considering model 
II for Pij and 3428 particles. 

• Figure 3: Growth velocity v along the y-direction proportional to | V^m | (curve 
A) and displacement | u \ (curve B). Curves C and D are for Poisson growth 
including screening and antiscreening, respectively. 

• Figure 4: Results for circular biharmonic growth using (a) model I with 2326 par- 
ticles and (b) model // with 2172 particles. Arrows indicate predicted transition 
points from Eq.(^. 

• Figure 5: Fractal nature of the simulated biharmonic patterns that are displayed 
in Figs. 4(a) and (b). Triangles are obtained using model / and squares using model 
//. Circles are for Laplacian growth giving a fractal dimension of ~ 1.7. 
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